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Abstract 

Simple lower and upper bounds on mean cycle time in stochas- 
tic acyclic fork-join queueing networks are derived using a (max, +)- 
algebra based representation of network dynamics. The behaviour of 
the bounds under various assumptions concerning the service times 
in the networks is discussed, and related numerical examples are pre- 
sented. 
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1 Introduction 

Fork-join networks introduced in [IJ [2], present a class of queueing system 
models which allow customers (jobs, tasks) to be split into several parts, and 
to be merged into one when they circulate through the system. The fork- 
join formalism proves to be useful in the description of dynamical processes 
in a variety of actual complex systems, including production processes in 
manufacturing, transmission of messages in communication networks, and 
parallel data processing in multi-processor computer systems. As a natural 
illustration of the fork and join operations, one can consider respectively 
splitting a message into packets in a communication network, each intended 
for transmitting via separate ways, and merging packets at a destination 
node of the network to restore the message. Further examples can be found 
in [J. 

The usual way to represent the dynamics of fork-join queueing networks 
relies on the implementation of recursive state equations of the Lindley type 
PQ. Since the recursive equations associated with the fork-join networks can 
be expressed only in terms of the operations of maximum and addition, there 
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is a possibility to represent the dynamics of the networks in terms of the 
(max,+)-algebra which is actually an algebraic system just supplied with the 
same two operations [31 El [5]. I* 1 fact) (max,+)-algebra models offer a more 
compact and unified way of describing network dynamics, and, moreover, 
lead to equations closely analogous to those in the conventional linear system 
theory [H El [8j [9]. In that case, the (max,+)-algebra approach gives 
one the chance to exploit results and numerical procedures available in the 
algebraic system theory and computational linear algebra. 

One of the problems of interest in the analysis of stochastic queueing 
networks is to evaluate the mean cycle time of a network. Both the mean 
cycle time and its inverse which can be regarded as a throughput present 
performance measures commonly used to describe efficiency of the network 
operation. 

It is frequently rather difficult to evaluate the mean cycle time exactly, 
even though the network under study is quite simple. To get information 
about the performance measure in this case, one can apply computer simula- 
tion to produce reasonable estimates. Another approach is to derive bounds 
on the mean cycle time. Specifically, a technique which allows one to estab- 
lish bounds based on results of the theory of large deviations as well as the 
Perron- Frobenius spectral theory has been introduced in |10j . 

In this paper we propose an approach to get bounds on the mean cycle 
time, which exploits the (max,+)-algebra representation of acyclic fork-join 
network dynamics derived in [HI [9] . This approach is essentially based on 
pure algebraic manipulations combined with application of bounds on ex- 
treme values, obtained in [TT| [T2]. 

The rest of the paper is organized as follows. Section [2] presents basic 
(max,+)-algebra definitions and related results which underlie the develop- 
ment of network models and their analysis in the subsequent sections. In 
Section [31 further algebraic results are included which provide a basis for 
derivation of bounds on the mean cycle time. 

A (max,+)-algebra representation of the fork-join network dynamics and 
related examples are given in Sectional Furthermore, Section [5] offers some 
monotonicity property for the networks, which is exploited in Section [6] 
to get algebraic bounds on the service cycle completion time. Stochastic 
extension of the network model is introduced in Section [71 The section 
concludes with a result which provides simple bounds on the network mean 
cycle time. Finally, Section [8J presents examples of calculating bounds and 
related discussion. 

2 Preliminary Algebraic Definitions and Results 

The (max, +) -algebra presents an idempotent commutative semiring (idem- 
potent semifield) which is defined as tR max = (JR> ©> <8>) with ]R = R U {e} , 
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e = — oo, and binary operations © and <8> denned as 



x © y = max(x, y), x(&y = x + y, for all x,y€.R. 

As it is easy to see, the operations © and © retain most of the prop- 
erties of the ordinary addition and multiplication, including associativity, 
commutativity, and distributivity of multiplication over addition. However, 
the operation © is idempotent; that is, for any x G JR, one has x © x = x. 

There are the null and identity elements in the algebra, namely e and 
, to satisfy the conditions x©e = e©x = x, and x®0 = 0©x = x, 
for any x G JR. The null element e and the operation © are related by the 
usual absorption rule involving x©e = e©x = £. 

Non-negative integer power of any x G R can be defined as x° = 0, and 
x©x <?_1 = x 9_1 ©x for q > 1. Clearly, the (max, +) -algebra power x q 



x H 



corresponds to qx in ordinary notations. We will use the power notations 
only in the (max, +) -algebra sense. 

The (max, +) -algebra of matrices is readily introduced in the regular 
way. Specifically, for any (n x n) -matrices X = (x^) and Y = {yij) , the 
entries of U = X © Y and V = X © Y are calculated as 



Uij — Xij © yij , and Vij =Q)x ik ®y kj . 

k=l 

As the null and identity elements, the matrices 

e ... e \ / £ \ 



I 



c ) 



0/ 



are respectively taken in the algebra. 

The matrix operations © and © possess monotonicity properties; that 
is, the matrix inequalities X < U and Y < V result in 



X © Y < U © V, 



X®Y < U®V 



for any matrices of appropriate size. 

Let X ^ £ be a square matrix. In the same way as in the conventional 
algebra, one can define X° = I , and X q = 10 X q ~ 1 = X q ~ x © X for any 
integer q > 1. However, idempotency leads, in particular, to the matrix 
identity 

(X © Y) q = X q ® X q ~ l © Y © • • • © Y q . 
As direct consequences of the above identity, one has 

(X © Y) q > x v © Y q - p , (I © X) q > (I © xy > x p , 

for all p = 0, 1, . . . , q. 
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For any matrix X, its norm is defined as 

ll^i^Tjl ~~ max Xi j . 

The matrix norm possesses the usual properties. Specifically, for any matrix 
X , it holds \\X\\ > e, and ||X|| = e if and only if X = £ . Furthermore, we 
have ||c(g>X|| = c(8) for any c G JR, as well as additive and multiplicative 
properties involving 

preY|| = ||x|| e ||y||, \\x®y\\ < \\x\\ \\y\\ 

for any two conforming matrices X and Y . Note that for any c > 0, we 
also have \\cX\\ = c||X||. 

Consider an (n x n) -matrix X with its entries Xij G JR. It can be treated 
as an adjacency matrix of an oriented graph with n nodes, provided each 
entry ^ e implies the existence of the arc in the graph, while 

= e does the lack of the arc. 

It is easy to verify that for any integer q > 1 , the matrix X q has its the 
entry 7^ e if and only if there exists a path from node i to node j in 
the graph, which consists of q arcs. Furthermore, if the graph associated 
with the matrix X is acyclic, we have X q = £ for all q > p, where p is the 
length of the longest path in the graph. Otherwise, provided that the graph 
is not acyclic, one can construct a path of any length, lying along circuits, 
and then it holds that X q 7^ £ for all q > . 

Consider the implicit equation in an unknown vector x = (xi, . . . ,x n ) T , 

x = U (g> x © v, (1) 

where U = (uij) and v = (v\, . . . , v n ) T are respectively given (n x n)- 
matrix and n-vector. Suppose that the entries of the matrix U and the 
vector v are either positive or equal to e. It is easy to verify (see, e.g. 
[31 Q2] that equation ([T]) has the unique bounded solution if and only if the 
graph associated with U is acyclic. Provided that the solution exists, it is 
given by 

x = (I(BU) P ® v, (2) 
where p is the length of the longest path in the graph. 

3 Further Algebraic Results 

We start with an obvious statement. 
Proposition 1. For any matrix X , it holds 

X < \\X\\ ® G, 

where G is the adjacency (e~0) -matrix of the graph associated with X. 
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Proposition 2. Suppose that matrices X\, . . . ,X& have a common associ- 
ated acyclic graph, p is the length of the longest path in the graph, and 

x = x^®.-.®x™\ 

where m\, . . . , mt are nonnegative integers. 

If it holds that m\ + • • • + m/, > p, then X = £ . 

Proof. It follows from Proposition [1] that 

X = X™ 1 © • • • © X™ k < ||Xi|| mi ® • • • <8> \\X k \\ mk © G mi+ - +mk , 

where G is the adjacency (e-O)-matrix of the common associated graph. 

Since the graph is acyclic, it holds that G q = £ for all q > p. Therefore, 
with q = mi + • • • + > p, we arrive at the inequality X < £ which leads 
us to the desired result. □ □ 

Lemma 1. Suppose that matrices X\, . . . ,Xf. have a common associated 
acyclic graph, and p is the length of the longest path in the graph. 
If \\Xi\\ > for all i = 1, . . . , k, then it holds 



(g)(/®Xir 



8=1 



< 



,i=l 



Hi 

k ! 



for any nonnegative integers mi, . . . , m^ . 
Proof. Consider the matrix 

k mi m k 

x = (gx/eLr^Qx; 1 ®-^^ 

i=l i'i=0 ifc=0 

mi m k 

= ©•••©^•••^ fc < © xi^...^x i 

ii=0 ifc=0 0<iiH hife<m 

where m = mi + • • • + m k ■ From Proposition [2] we may replace m with p 
in the last term to get 

x< xi 1 ®--®x*. 

0<iiH Hk<P 

Proceeding to the norm, with its additive and multiplicative properties, 
we arrive at the inequality 

11*11 < prill* 1 ®-"®!!**!!*. 

0<iiH Mfc<p 

Since for all i = 1, . . . , k, it holds < ||*|| < ||*i|| © • • • © ||*fc|| j we finally 
have p 

||*||<©(||*i||ffi--ffi||* fc ||) p = (©11*11) • 

i=0 \i=0 J 

□ □ 
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4 An Algebraic Model of Queueing Networks 

We consider a network with n single-server nodes and customers of a sin- 
gle class. The topology of the network is described by an oriented acyclic 
graph Q = (N, A) , where the set N = {1, . . . , n} represents the nodes, and 
A = {(hj)} C N x N does the arcs determining the transition routes of 
customers. 

For every node i G N , we denote the sets of its immediate predecessors 
and successors respectively as P(i) = {j\ G A} and S(i) = {j\ G 
A}. In specific cases, there may be one of the conditions P(i) = and 
S(i) = encountered. Each node i with P(i) = is assumed to represent 
an infinite external arrival stream of customers; provided that S(i) = 0, 
it is considered as an output node intended to release customers from the 
network. 

Each node i G N includes a server and its buffer with infinite capacity, 
which together present a single-server queue operating under the first-come, 
first-served (FCFS) discipline. At the initial time, the server at each node i 
is assumed to be free of customers, whereas in its buffer, there may be r^, 
< Ti < oo, customers waiting for service. The value v% = oo is set for 
every node i with P(i) = 0, which represents an external arrival stream of 
customers. 

For the queue at node i, we denote the kth. arrival and departure epochs 
respectively as Ui(k) and Xi(k) . Furthermore, the service time of the kth 
customer at server i is indicated by . We assume that > are given 
parameters for all i = 1, . . . ,n, and k = 1,2,..., while Ui(k) and Xi(k) are 
considered as unknown state variables. With the condition that the network 
starts operating at time zero, it is convenient to set Xj(0) = 0, and Xi{k) = e 
for all k < , i = 1 , . . . , n . 

It is easy to set up an equation which relates the system state variables. 
In fact, the dynamics of any single-server node i with an infinite buffer, 
operating on the FCFS basis, is described as 

Xi{k) = r ik (8) Ui(k) T ik ®Xi(k - 1). (3) 



With the vector-matrix notations 





( Mk) \ 




( Mk) \ 




( T\k 


£ 


u(k) = 




, x(k) = 




, Tk = 








\ u n (k) ) 




\ x n (k) J 









we may rewrite equation ([3]) in a vector form, as 

x{k) =T k ® u{k)@T k ® x(k- 1). (4) 
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4.1 Fork- Join Queueing Networks 

In fork-join networks, in addition to the usual service procedure, special join 
and fork operations are performed in its nodes, respectively before and after 
service. The join operation is actually thought to cause each customer which 
comes into node i, not to enter the buffer at the server but to wait until 
at least one customer from every node j G P(z) arrives. As soon as these 
customers arrive, they, taken one from each preceding node, are united into 
one customer which then enters the buffer to become a new member of the 
queue. 

The fork operation at node i is initiated every time the service of a cus- 
tomer is completed; it consists in giving rise to several new customers instead 
of the original one. As many new customers appear in node % as there are 
succeeding nodes included in the set S(i). These customers simultaneously 
depart the node, each being passed to separate node j G S(z). We assume 
that the execution of fork-join operations when appropriate customers are 
available, as well as the transition of customers within and between nodes 
require no time. 

As it immediately follows from the above description of the fork-join 
operations, the kth arrival epoch into the queue at node i is represented as 

Ui{k) = < jeP(i) (5) 
I e, if P(i) = 0. 

In order to get this equation in a vector form, we first define the number 
M = max{rj| T{ < oo, i = 1, . . . , n} . Now we may rewrite ([5|) as 

M n 

Mb) = ®®97i®Xj{k-m), 

m=0 j=l 



where the numbers g™ are determined by the condition 

(6) 



0, if i G P(j) and m = r,-, 



e, otherwise. 



Let us introduce the matrices G m = yg^j for each m = 0, 1, . . . , M . In 
fact, G m presents an adjacency matrix of the partial graph Q m = (N, A TO ) 
with A m = G A; Tj = m} . Since the graph of the entire 

network is acyclic, all its partial graphs Q m , m = 0, 1, . . . ,M, possess the 
same property. 

With these matrices, equation (|5|) may be written in the vector form 

M 

u(k) = 0G^(g) x(k-m), (7) 

m=0 
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where denotes the transpose of the matrix G m . 

By combining equations @) and ©, we arrive at the equation 

x(k) = T k ®Gl® x(k) T k <8> x{k - 1) 

M 

®T k ®(^)G T l ®x(k-m). (8) 

m=l 

Clearly, it is actually an implicit equation in x(k), which has the form of 
(UJ), with U = Tk®GQ . Taking into account that the matrix Tk is diagonal, 
one can prove the following statement (see also [SJE]). 

Theorem 2. Suppose that in the fork- join network model, the graph Qq 
associated with the matrix Gq is acyclic. Then equation JS^ can be solved 
to produce the explicit dynamic state equation 

M 

x ( k ) = A m(k)<S> x(k-m), (9) 

m=l 

with the state transition matrices 

Ai(fc) = {I®T k ®GlY®Tk® (/eGf), (10) 
A m (k) = {I®T k ®Gl) p ®T k ®G T m , m = 2,...,M, (11) 

where p is the length of the longest path in Qq . 
4.2 Examples of Network Models 

An example of an acyclic fork-join network with n = 5 is shown in Fig. [TJ 
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Figure 1: An acyclic fork-join network. 
Since for the network M = 1 , we have from (0) 
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Taking into account that for the graph Qq , the length of its longest path 
p = 2, we arrive at the dynamic equation 

x(k) = A(k) ®x{k- 1), 

with the state transition matrix calculated from (|1U|) as 

A(k) = (I@T k ® Glf ®T k ®{I@ Gl) 



I n fe e e e e \ 

e r 2 k e e e 

Tlk®T?, k £ T 3 fc £ £ 

r 4fc r 4fc £ r 4 fc E 



\ (Tlfc<g>T 3fc eT4fc)<g>T5fc T 4k ®T 5k T 3k ®T 5k T 4k ®T 5k T 5k / 

Note that open tandem queueing systems (see Fig. [2]) can be consid- 
ered as trivial networks in which no fork and join operations are actually 
performed. 

12 n 

□o — do — do — - 

n = oo T2 = r n = 



Figure 2: Open tandem queues. 



For the system in Fig. [2j we have M = , and p = n — 1 . Its related 
state transition matrix A(k) has the entries [SJ [7] 



atj(k) 



Tjk®T j+lk ®---®T ik , if i>j, 

e, otherwise. 



5 A Monotonicity Property 

In this section, a property of monotonicity is established which shows how 
the system state vector x (k) may vary with the initial numbers of customers 
T{ . It is actually proven that the entries of x (k) for all k = 1, 2, . . . , do not 
decrease when the numbers rj with 0<rj<oo, i = 1,... ,n, are reduced 
to zero. 

As it is easy to see, the change in the initial numbers of customers results 
only in modifications to partial graphs Q m and so to their adjacency matrices 
G m . Specifically, reducing these numbers to zero leads us to new matrices 
G = Go © Gi ■ ■ ■ G M , and G m = £ for all m = 1, . . . , M. 

We start with a lemma which shows that replacing the numbers = 1 
with ri = does not decrease the entries of the matrix A\{}z) defined by 

(Hod. 
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Lemma 3. For all k = 1, 2, . . . , it holds 

A x {k) < A(k) 

with A(k) = (I © Tk © G^)i ® Tk , where Gq = Go © G\ , and q is the length 
of the longest path in the graph associated with the matrix Go . 

Proof. Consider the matrix A\{k) and represent it in the form 

Ai(fc) = ((I © T k © Gq) p ® T k ) © ((/ © T k © GTf © % © Gf ), 

where p is the length of the longest path in the graph associated with Go • 
As one can see, to prove the lemma, it will suffice to verify both inequal- 
ities 

A(k) > {I®%®Gl) p ®T k , (12) 
A(k) > (I®T k ®G r ) p ®T k ®Gj. (13) 

Let us write the obvious representation 



{I@T k ® GqY = 0(7 © T k ® GlY ® {T k © Gj 



yT\q-i 
i=0 

Since q > p, we get from the representation 

{I@T k ®Gl) q > {I®T k ®Glf = {{I®T k ®Gl)®T k ®Gly> > (I®T k ®G%) p . 

It remains to multiply both sides of the above inequality by Tk on the right 
so as to arrive at (|12p . 

To verify (fTB"|) . let us first assume that q > p. In this case, we obtain 

(I © % © GlY >(I®T k ® GlY +l 

= ((I © T k © Gft) © % © Gf > (/ © T fc © Glf © T fc © Gf. 

Suppose now that q = p. Then it is necessary that G\ © Gq = £ . If this 
were not the case, there would be a path in the graph associated with the 
matrix Go = Go © G\ , which has its length greater than p , and we would 
have q > p. 

Clearly, the condition G\ © Gq = £ results in {T k © Gq ) p © Tk © Gj = £ , 
and thus we get 

(/ © % © Gq ) 9 = ((I © T k ® Go") © T fc ® Gf ) p 

> (l0 7i®Gg , ) p - 1 ®'ft®Gf = (/©r fc ©G r ) p ©r fc ©Gf. 

Since it holds {I@T k ® G%) p ®T k >{I@T k ® G%) p , one can conclude that 
inequality (fT3|) is also valid. □ □ 
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Theorem 4. In the acyclic fork-join queueing network model iTP Hiij) . re- 
ducing the initial numbers of customers from any finite values to zero does 
not decrease the entries of the system state vector x(k) for all k = 1,2, ... . 

Proof. Let x(k) be determined by P HTTj) . Suppose that the vector x(k) 
satisfies the dynamic equation 

x(k) = A(k) (8) x(k - 1) 

with 

(M \ q 
I®T k ®($G T m \ ®T k = {I@T k ®G T y®T k , 

where q is the length of the longest path in the graph associated with the 
matrix G = G Gi • • • G m . 

Now we have to show that for all k = 1,2, ... , it holds 

x(k) < x(k). 

Since x(k\) < x(k2) for any k\ < k2, we have from 

M / M \ 

x(k) = 0^ m (A;)0 x(k-m) < I 0i m (fc) J x(k-l). 

m=l \m=l J 

Consider the matrix 

M / M \ 

M(k) = ®A m {k) = {I®T k ®G%Y®%® I iffi0Gl • 

m=l \ m=l / 

By applying Lemma [31 we have 

/ M \ q 

M{k)< /er fe ®Gj®0GU ®r fe = 2(fe). 

V m=l / 

Starting with the condition x (0) = £ (0) , we successively verify that the 
relations 

x(k) < Ai(k) ® x(k-l)< A(k) ®x(k-l)< A(k) ®x(k-l) = x(k) 
are valid for each k = 1, 2, □ □ 
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6 Bounds on the Service Cycle Completion Time 

We consider the evolution of the system as a sequence of service cycles: the 
1st cycle starts at the initial time, and it is terminated as soon as all the 
servers in the network complete their 1st service, the 2nd cycle is terminated 
as soon as the servers complete their 2nd service, and so on. Clearly, the 
completion time of the kth cycle can be represented as 

maxxi(k) = \\ x (k)\\. 

i 

The next lemma provides simple lower and upper bounds for the kth 
cycle completion time. 

Lemma 5. For all k = 1, 2, . . . , it holds 
k 



k / k \ 

< Era +p ©ii*ii 

i=l \i=l / 



< \\X( 

i=l 

Proof. To prove the left inequality first note that 

Ai(k) = (1 © T k © GlY ® T k Gj) > T k . 

With this condition, we have from ([9]) 

M 

x(k) = (J) A m (k) (g> x(k-m)> A x {k) <g) x(k - 1) > Tk ® x{k - 1). 

m=l 

Now we can write 

x(k) >T k ®x{k-l) >Tk®T k -x®x{k-2) >•• > % ® ■ ■ ■ <g> Ti ® x(0), 

where x(0) = 0. Taking the norm, and considering that 71, i = 1, ...,k, 
present diagonal matrices, we get 

\\x(k)\\ > \\T k ® •••®71|| = \\Ti + ■■■ + T k \\. 

To obtain an upper bound, let us replace the general system p HTTj) with 
that governed by the equation 

x(k) = A(k) ®x(k — l) (14) 

with A(k) = (I ®T k ®G T ) q ®T k , where G = G © Gi ® ■ ■ ■ G m , and q is 
the length of the longest path in the graph associated with G . As it follows 

from Theorem [31 one has x(k) < x(k) for all k = 1, 2, 

Let us denote Ak = A(k) ® • • • (g> A(l). With the condition x(0) = 
x(0) = , we get from ([1] 



\x(k)\\ = \\A(k) ® ■ ■ ■ ® A(l)\\ = \\A h 
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With Proposition [2] we have 



A k = (g)(J © T k -i + i ® G T ) q © %- i+ i < (g) ||75|| © (g)(I © r fc - i+ i © G T ) q . 

i=l i=l i=l 

Proceeding to the norm and using Lemma [TJ we arrive at the inequality 

k / k \ q 



\A k \\ <Q<)\\Ti\ 
1=1 



' k \ q k / k \ 

era =Era+W©ra ■ 

vi=l / i=l \i=l / 



which provides us with the desired result. □ □ 



7 Stochastic Extension of the Network Model 

Suppose that for each node i = 1, . . . , n, the service times Tji,Tj2, • • • , form 
a sequence of independent and identically distributed (i.i.d.) non-negative 

random variables with E[th~] < oo and D[rjfc] < oo for all k = 1, 2, 

As a performance measure of the stochastic network model, we consider 
the mean cycle time which is defined as 

7 = lim Ux(k)\\ (15) 
fc-s>oo k 

provided that the above limit exists. Another performance measure of in- 
terest is the throughput defined as ir = I/7 . 

Since it is frequently rather difficult to evaluate the mean cycle time 
exactly, even though the network under study is quite simple, one can try 
to derive bounds on 7. In this section, we show how these bounds may be 
obtained based on (max,+)-algebra representation of the network dynamics. 

We start with some preliminary results which include properties of the 
expectation operator, formulated in terms of (max,+)-algebra operations. 



7.1 Some Properties of Expectation 

Let £1,... be random variables taking their values in R, and such that 
their expected values E[£j] , i = 1, . . . , k, exist. 

First note that ordinary properties of expectation leads us to the obvious 
relations 



©6 

,i=i 



<®E&], 
i=l 



and 



i=l 



i=i 



Furthermore, the next statement is valid. 
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Lemma 6. It holds 



E 



06 

.i=i 



i=l 



Proof. The statement of the lemma for k = 2 follows immediately from the 
identity 

x ® y = -(x + y + | a — y|), forallx,y€lR 

and ordinary properties of expectation. It remains to extend the statement 
to the case of arbitrary k by induction. □ □ 

The next result [HI [12] provides an upper bound for the expected value 
of the maximum of i.i.d. random variables. 

Lemma 7. Let £i, . . . be i.i.d. random variables with E[£i] < oo and 
D[£i] < oo. Then it holds 



Consider a random matrix X with its entries taking values in JR. 
We denote by E\X] the matrix obtained from X by replacing each entry 
by its expected value E[xjj] . 

Lemma 8. It holds 

E||X||>||E[X]||. 
Proof. It follows from Lemma [6] that 



EIIXII = E 



>0E[x, 



|E[X]||. 



□ 



□ 



7.2 Existence of the Mean Cycle Time 

In the analysis of the mean cycle time, one first has to convince himself that 
the limit at (|15p exists. As a standard tool to verify the existence of the 
above limit, the next theorem proposed in [14J is normally applied. One can 
find examples of the implementation of the theorem in the (max, +) -algebra 
framework in [1CH 0] . 

Theorem 9. Let {£ik\ I, k = 0, 1, . . . ; I < k} be a family of random variables 
which satisfy the following properties: 

Subadditivity: < £j m + £ m k for all I < m < k ; 
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Stationarity: both families {£i+ik+i\ I < k} and {£ik\ I < k} have the 
same joint distributions; 

Boundedness: for all k = 1,2,..., there exists E[£ofc] > —ck for some 
finite number c. 

Then there exists a constant 7, such that it holds 

1. lim £,ok/k = 7 with probability 1, 

2. lim E[£o*]/fc = 7- 

k— >oo 

For simplicity, we examine the existence of the mean cycle time for a 
network with the maximum of the initial numbers of customers in nodes 
M < 1. As it follows from representation (f9 Hlip . the dynamics of the 
system may be described by the equation 

x(k) = A(k) (g> x(k - 1) 

with the matrix A(k) = Ai(k) determined by (|10p . Clearly, in the case 
of M > 1, a similar representation can be easily obtained by going to an 
extended model with a new state vector which combines several consecutive 
state vectors of the original system. 

To prove the existence of the mean cycle time, first note that with k = 
1,2, ... , are i.i.d. random variables for each i = 1, . . . ,n, and consequently, 
7fe are i.i.d. random matrices, whereas [|7fc|| present i.i.d. random variables 
with E||7fc|| < 00 and D[|7fc[| < 00 for all A; = 1,2,.... 

Furthermore, since the matrix A(k) depends only on 7fc, the matrices 
A(1),A(2), . . . , also present i.i.d. random matrices. It is easy to verify that 
< E||A(Jfe)|| < 00 for all k = 1,2,.... 

In order to apply Theorem [9] to stochastic system Q with transition 
matrix (|10|) , one can define the family of random variables | I < k} with 

b k = \\A(k)®---®A(i + i)\\. 

Since A(i) , i = 1,2,..., present i.i.d. random matrices, the family 
{£lk\ I < k} satisfies the stationarity condition of Theorem [9l Furthermore, 
the multiplicative property of the norm endows the family with subaddi- 
tivity. The boundedness condition can be readily verified based on the 
condition that < E[tj/%] < 00 for all i = 1, . . . ,n, and k = 1,2, . . . . 

7.3 Calculating Bounds on the Mean Cycle Time 

Now we are in a position to present our main result which offers bounds on 
the mean cycle time. 

Theorem 10. In the stochastic dynamical system (OJ) the mean cycle time 
7 satisfies the double inequality 

||E[Ti]||<7<E||Ti||. (16) 
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Proof. Since Theorem [9] hods true, we may write 

1 

7 = lim — E||aj(fc)||. 

k— too k 

Let us first prove the left inequality in (|16|) . From Lemmas [5] and [8l we 

have 



-Ellcc 
k 



> 



i 



i=l 



|E[Ti]||, 



independently of k . 

With the upper bound offered by Lemma [5l we get 



^E||a;(A!)|| <E||7I|| + |E 



i=i 



From Lemma [?1 the second term on the right-hand side may be replaced by 
that of the form 



?(F|7i| 



fc- 1 

y/2k - 1 



VDHTill , 



which tends to as k — > oo . 



□ 



□ 



8 Discussion and Examples 

Now we discuss the behaviour of the bounds (|16|) under various assumptions 
concerning the service times in the network. First note that the derivation 
of the bounds does not require the fcth service times to be independent 
for all i = 1, . . . , n. As it is easy to see, if = for all i, we have 
||E[7~]i|| = E||71||, and so the lower and upper bound coincide. 

To show how the bounds vary with strengthening the dependency, we 
consider the network with n = 5 nodes, depicted in Fig. [TJ Let th = 
Y^=i a ijCjii where j = 1,...,5, are i.i.d. random variables with the 
exponential distribution of mean 1 , and 

_ J a, if % = j, 

a ^-\ 1(1 -a), ifi^j, 

where a is a number such that 1 < a < 1/5. 

It is evident that for a = 1 , one has t& = £n , and then th, i = 1, . . . , 5 , 
present independent random variables. As a decreases, the service times th 
become dependent, and with a = 1/5, we will have Tn = (£n + • • • + ^5i)/5 
for all i= 1, . . . , 5. 

Table [1] presents estimates of the mean cycle time 7 obtained via simula- 
tion after performing 100000 service cycles, together with the corresponding 
lower and upper bounds calculated from (fTE]). 
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a 


l|E[Ti]|| 


7 


E||7i|| 


1 


1.0 


1.005718 


2.283333 


1/2 


1.0 


1.002080 


1.481250 


1/3 


1.0 


1.000871 


1.213889 


1/4 


1.0 


1.000279 


1.080208 


1/5 


1.0 


1.000000 


1.000000 



Table 1: Numerical results for a network with dependent service times. 



Let us now consider the network in Fig. [T] under the assumption that the 
service times th are independent exponentially distributed random vari- 
ables. We suppose that E[rji] = 1 for all i except for one, say i = 4, with 
E[t4i] essentially greater than 1. One can see that the difference between 
the upper and lower bounds will decrease as the value of E[t4i] increases. 
Table [2] shows how the bounds vary with different values of E[t4i] . 



E[T4l] 


l|E[7i]|| 


7 


E||7iJ| 


1.0 


1.0 


1.005718 


2.283333 


2.0 


2.0 


2.004857 


2.896032 


3.0 


3.0 


3.004242 


3.685531 


4.0 


4.0 


4.003627 


4.554525 


5.0 


5.0 


5.003013 


5.465368 


6.0 


6.0 


6.002398 


6.400835 


7.0 


7.0 


7.001783 


7.351985 


8.0 


8.0 


8.001168 


8.313731 


9.0 


9.0 


9.000553 


9.282968 


10.0 


10.0 


10.000008 


10.257692 



Table 2: Results for a network with a dominating service time. 

Let us discuss the effect of decreasing the variance D[rji] on the bounds 
on 7 . Note that if th were degenerate random variables with zero variance, 
the lower and upper bounds in (fT6"j) would coincide. One can therefore 
expect that with decreasing the variance of th , the accuracy of the bounds 
increases. 

As an illustration, consider a tandem queueing system (see Fig. [2]) with 
n = 5 nodes. Suppose that m = £a/r, where , i = 1, . . . , 5, are i.i.d. 
random variables which have the Erlang distribution with the probability 
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density function 



/rW ~\ 0, if t < 0. 

Clearly, E[th] = 1 and E)[th] = 1/r. Related numerical results including 
estimates 7 evaluated by simulating 100000 cycles are shown in Table [3j 



r 


l|E[7I]|| 


7 


E||Ti|| 


1 


1.0 


1.042476 


2.928968 


2 


1.0 


1.026260 


2.311479 


3 


1.0 


1.019503 


2.045538 


4 


1.0 


1.015637 


1.890824 


5 


1.0 


1.013110 


1.787242 


6 


1.0 


1.010864 


1.711943 


7 


1.0 


1.009920 


1.654154 


8 


1.0 


1.008409 


1.608064 


9 


1.0 


1.007726 


1.570232 


10 


1.0 


1.006657 


1.538479 



Table 3: Results for tandem queues at changing variance. 
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